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Abstract 

In spallation neutron sources, liquid mercury is the subject of big thermal and 
pressure shocks, upon adsorbing the proton beam. These changes can cause unstable 
bubbles in the liquid, which can damage the structural material. While there are 
methods to deal with the pressure shock, the local temperature shock cannot be 
avoided. In our paper we calculated the work of the critical cluster formation (i.e. 
for mercury micro-bubbles) together with the rate of their formation (nucleation 
rate). It is shown that the homogeneous nucleation rates are very low even after 
adsorbing several proton pulses, therefore the probability of temperature induced 
homogeneous bubble nucleation is negligible. 

1 Introduction 

Irradiating liquid metal (usually mercury) with proton beams is up to now the best 
method to produce high-intensity, multi-purpose neutron beams. This method has been 
used in various existing facilities and it is planned to be used in the European Spallation 
Source, too. Unfortunately upon adsorbing the high-intensity proton beam in the liquid, 
the neutrons are not the only ones emitted; an unavoidable heat and pressure wave will 
be emitted simulataneously from the adsorption region. The increase of the temperature 
and (in the negative period of the pressure wave) the decrease of the pressure can cause 
cavitation in the liquid. The metal vapor bubbles then will flow with the liquid and upon 
reaching high pressure and low temperature regions, they will collapse, causing severe 
damage in nearby solid structures. This phenomenon is known as cavitation erosion 
and one of the main factors which (due to pitting and weight loss) shorten the lifetime 
of structural materials significantly. Therefore to avoid cavitation is one of the main 
challenges of the design of the spallation source target [H |2l [3]. 

It should be mentioned here, the along the methods to minimize cavitation itself, there 
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are two other ways to minimize the damage. One of them are the various ways of surface 
treatments (plasma nitriding, plasma carbonizing, etc.), which makes the surface more 
resistant to the damaging pressure wave emitted by the collapsing bubble [U [5]. The 
other one is the addition of helium micro-bubbles, which is a proven way to soften up and 
to reduce the damage by absorbing the expansion of liquid mercury and mitigating the 
pressure waves [H El [6l [71 18] . Considering this method as a successful one to deal with 
the pressure-drop induced cavitation, in our paper we focused mainly on the temperature 
increase induced cavitation and allowed only small pressure changes to occur (down to -5 
bar). 

Our main aim is to determine whether the conditions (temperature and pressure 
changes) are able to cause cavitation or not. We approached the problem in three steps. 
In the first step (Section [2]), we calculated the phase equilibrium, stability limit and var- 
ious other properties of mercury by using the equation of state proposed by Morita et al. 
[HI [ini [m 112] • In the next step (Section [3|), we made some estimation for the magnitude 
of pressure and temperature changes by using single and repeated proton pulses on mer- 
cury. In the final step (Section H]), we calculated the work of critical bubble formation in 
mercury as well as the rate of homogeneous nucleation in the pressure-temperature range 
defined according to the results of the previous section. The paper is completed by a 
summary and discussion (Section [5|). 

2 Model system 

2.1 Location of binodal and spinodal curves 

For the description of mercury (Hg) in both the liquid and gas phases, we will apply a 
slightly modified thermal equation of state as compared to the expression proposed by 
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Morita et al. (see [HI [TO], and, in particular, Eq. (15) in [H]). It reads 

p = ^^ ^ (1) 

^ K{T){v-b) v{v + c)' 

a(T) = a, (^-j at T < T, , (2) 

where R = 8.314 J • moP^ is the universal gas constant, p is pressure, v is molar volume, 
T is temperature, Oc, b, c and n are the model parameters specific for the substance 
considered, Tc is critical temperature. The correction coefficient K{T) is dependent on 
temperature only, it was introduced in the repulsive term instead of the parameter Xd, 
which is a function of T and p (see [12], in such case Eq. ([1]) becomes an equation for 
definition of p{v,T), and has no analytical solution). 
We employ further dimensionless variables 

n = ^, u; = -, e = ^, (3) 

Pc Vc Tc 

where Vc is the molar volume, pc the pressure both at the critical point with the critical 
temperature, Tc. These parameters can be determined from Eq. ([1]) in the common way 



via 

/ dp\ f d'^p\ 



at T = Tc. (4) 



The equation of state in reduced variables is given by 



Here 

^J0\ = 

RTc 



xM = (6) 



is the reduced critical compressibility, and 

ir(g) = 1.106697 - 0.106697 -exp I V ^ 1 , (7) 



PcV^ 



exp j 


f e-1 \ 


^0.17026/ 


b 






e = - 


Vc 


Vc 



«W = ^ = «^", /9 = -, e = -- (8) 



According to we have then 

a = 2.5272, /3 = 0.3952 , ^ = -0.16567 , n = -0.0284127 . (9) 
From Eqs. ([I]) and @ we get [11] 

We = 1.797 • 10"^ mVkg , Pc = 5566 kg/ia , (10) 
= 158 • 10*^ Pa , Te = 1762 K . 

The location of the classical spinodal curve can be found via the determination of the 
extrema of the thermal equation of state, 11(6', u) (Eq. (Q) considering the temperature 9 
as constant. By taking the derivative of n(^, u) with respect to oj, we obtain from Eq. ([5j) 
the result 

d aimu^ 9 _ 

^ ' ^ " WTeF " xc{u;-py - ° • ^^^^ 

For 9 < 1, this equation has two positive solutions co^^^^^ and w^p^^*-* for co corresponding 
to the specific volumes of the both macrophases at the spinodal curves (or at the limits 
of metastability). 

Similarly, the binodal curves give for 9 < 1 the values of the specific volumes of 
the liquid and the gas phases coexisting in thermal equilibrium at a planar interface. 
From the left branch of the binodal curve, we get the specific volume of the liquid phase 
{uji^'^'^\9) = uj^'^^^\9)), from the right branch of the binodal curve, we obtain the specific 
volume of the gas {u^^^\9) = ujI^^^^^\9)) . For 9 = 1, both solutions coincide in the critical 
point {ool'^'^^ = oojj^^'^ = ujc = 1), again. Consequently, in order to determine the specific 
volumes of the liquid and the gas at some given temperature in the range ^ < 1, we have 
to specify the location of the binodal curve. 

The location of the binodal curve may be determined from the necessary thermody- 
namic equilibrium conditions (for planar interfaces) - equality of pressure and chemical 
potentials - via the solution of the set of equations 

Ui{ui, 9) = Yigiuog, 9) , pi{uou 9) = pg{ujg, 9) . (12) 



Here by fi the chemical potential of the atoms or molecules in the liquid (/) and the gas 
(g) are denoted. Having at our disposal already the equation for the reduced pressure 
(c.f. Eq. ([5])), we have now to determine in addition the chemical potential in dependence 
on pressure and temperature. This task will be performed in the next section. 

Isotherms for mercury according to Eq. for different values of the reduced temper- 
ature 9 = 0.4, 0.65, 0.8, 0.891 and 0.92 are shown in Fig. [H dashed and dashed-dotted 
curves present binodal and spinodal, correspondingly. One can see, that there are two 
classes of isotherms: for the first one {6 > 0^) p > 0, and for the second class {9 < 9^) 
pressure may be both positive and negative. The parameter 9s is determined via the 
equation 



for mercury 9s ~ 0.891 and T, = Tc9s ~ 1570 K. A comparison of experimental data 
[T3| [H] for the vapor-liquid coexistence properties of mercury with results obtained in 
this work are shown in Fig. |2] and Fig. [3] for (T, p) and (p, T~^)-variables, correspondingly. 

2.2 Determination of the chemical potential and the interfacial 
tension 

For processes at constant temperature, the change of the Helmholtz free energy, F, may 
be expressed as 



Here V is the volume of the system and n the number of moles in it. For a given fixed 
mole number, n, of the substance (n = constant), we have, in particular. 



IiiM9s),9s) = Q, 



(13) 



dF = —pdV + ndn . 



(14) 



dip, 



—pdv , 



F 



V 

n 



(15) 



'n 



n 



or, in reduced variables. 




(16) 
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Employing in the integration of Eq. f|T6|) the equation of the state, Eq. ([5]), we obtain 



,PcVc 



aid) f C\ ^ , 
-^In 1 + ^ +^lnu;-/3 



(17) 



Alternatively, the change of the Helmholtz free energy - provided the volume V is 
fixed - is given at constant temperature by 



dF = fidn 



(18) 



From Eq. ( !T8|) . we arrive at 



d(pi 



~dv , 

V 



(19) 



On the other side, the functions ip^ and ip^ are connected by 



F = ipnU = ip^V 



Vv = 

V 



(20) 



With Eq. ( IT7|) . we have then 



Pc 

'^v = — 



"^'^nfl + iU4^1n(.-/5) 



(21) 



With Eqs. (fT9|) and ( 1211) . the expression for the chemical potential of a HLM can be 
obtained then via 

This relation yields 



PcVc 



a{9) 



Quo 



+ ^lnfl + iU^ln(.-/5) 



(23) 



In addition to the bulk properties of the system under consideration, we have to know 
the value a of the surface tension for a coexistence of both phases at planar interfaces in 
dependence on the parameters describing the state of both phases. We choose here this 
dependence in the form [151 IISI El HE] 



1 1 



6 = 2.5 



(24) 



where 



e{9) = A 



1 1 

(irfty ~ (right) 



n—5 



(25) 



and A and n are constant parameters. Comparison of Eqs. flM|) and fl25l) with experimental 
data [19] 

(t(T) = 0.5446544 - 0.000204917 ■ T (26) 

(vahd in this form only for temperatures far below the critical temperature; here the 
temperature is given in Kelvin and the surface tension in J/m^) at cui = w^'*^^*^ and Ug = 



(right) 



yields 



A = 0.033253 J/in 



n = 3 . 



(27) 



In Fig. m dependence of the surface tension on temperature is shown, solid curve presents 
Eq. at uJi = uj^'^^^\ ujg = u;^"^'^*'', and dashed curve - Eq. 



3 Determination of the pressure and temperature 
change after proton adsorption 

For the determination of the pressure and temperature change, a "one dimensional six- 
equation two-fluid model" was used, which is capable to describe transients like pressure 
waves, quick evaporation or condensation which is proportional to cavitation caused by 
energetic proton interaction in mercury target |2U]. The method was developed to describe 
the sudden and drastic steam condensation, called water hammer [2T| [22]. 

The model contains six first-order partial-differential equations which describe one- 
dimensional surface-averaged mass, momentum and energy conservation laws for both 
phases. A special numerical procedure ensures that shock-waves can be described without 
any numerical dispersion. With two major modifications this model can be applied to 
investigate the thermo-hydraulic properties of the planned mercury target in the European 



8 



Spallation Source (ESS). These modifications are the following: the equation of state 
namely the density and the internal energy of both mercury phases should be known in 
a broad range of pressure (1 Pa to 100 MPa) and temperature (273 K to 1000 K). As 
a second point the interaction of the high energy proton beam with mercury has to be 
included. This is a much simpler task because we may consider that about 50 % of the 
300 kJ/pulse beam energy is absorbed as a 2ms long heat shock square pulse, giving a new 
source terms in the energy equation of the liquid phase. The ESS mercury target station 
is modeled as a 18 meter long closed loop which is in three dimension the pipe diameter 15 
cm. We consider that 150 kJ heat is absorbed in a 10 cm long pipe, this is approximately 
the width of the proton pulse. Calculation shows that such a single pulse heats up the 
mercury with about 40-44 K, assuming that the initial temperature was between 293-373 
K (i.e. within the normal working range of the spallation source). In the calculations, low 
velocities 0.5-4 m/s, low initial pressure 1-4 bar and low initial temperature (below 374 
K) were assumed. To our knowledge the existing Japanese Spallation Neutron Source Hg 
loop is about 15 m long, with a diameter of 15 cm, the flow velocity of Hg is 0.7 m/s and 
the pressure is approximately equal to 1 bar. 

Concerning the pressure change, the model is able to estimate the positive part, but 
at the negative region (where most of the low temperature cavitation is expected to 
happen [23j), a stability problem aroused. Therefore we focused our calculation to the 
heat shock and, at present, neglected the pressure change. Preliminary calculation yielded 
a few bar changes [20], in agreement with the results of Ida [211 125], therefore the latter 
calculations were performed in the -5 to 10 bar range. We should mention here, that other 
models predicted much larger pressure changes (even hundreds of bars) EZ] both in 
the positive and negative pressure region. 

Also the effects of repeated pulses were checked. The calculations were performed 
with a 2 ms square pulse train where the delay time was 20 ms which is similar to a 16Hz 



9 



repetition rate. We started with a flow system with initial p=4 bar, Tmiuai = 353 K and 
initial flow velocity v = 4 m/s. We found that the temperature jumps are more or less 
additive which means that after the beginning of the third pulse the temperature was 
about 430 K. 

4 Determination of the work of critical cluster for- 
mation 

Let us assume, now, that the system is brought suddenly into a metastable state located 
between binodal curve and spinodal curve at the liquid branch of the equation of state. 
Then, by nucleation and growth processes, bubbles may appear spontaneously in the 
liquid and a phase separation takes place [28]. Based on the relations outlined above, we 
will determine now the parameters of the critical clusters governing bubble nucleation in 
dependence on the state parameters, pressure and temperature. 

We start with the general expression for the change of the thermodynamic potential 



Here the subscript a specifies the parameters of the cluster (bubble) phase while (3 refers 
to the ambient liquid phase. This relation holds generally provided - as we assume - the 
state of the ambient liquid phase remains unchanged by the formation of one bubble. For 
a one-component system (as discussed here), this expression is reduced to 



As the independent variables, we select the size of the bubble, r and the molar volume of 
the gas phase in the bubble. Similarly to [151 ISH EQ], we arrive then at 



(28) 



j 



AG = aA+ {p- Pa)Va + n„ [fi^ - /i/3] . 



(29) 




(30) 



10 



where the following notations have been introduced: 

/ K, 9) = nK, 9) - Uiu, 9) + ^ f ,M-,M \ ^3^^ 

(jjg Y Pc^c ) 

9^^, ^1= ^-^^^MO? , (32) 
ill 3 pikBTc9 

r = ^ , K = -e{9) . (33) 

Ra Pc 



The dependence of the scaling parameters Qi and R^ on the reduced temperature is shown 
in Fig. El 

The Gibbs free energy surface for the metastable initial state has typical saddle shape 

at the critical point (see Fig. ^ 9 = 0.92, ui = 0.65). The critical point position is 

determined by the set of equations 

dAgjr, ujg,uji,9) ^ ^ dAg{r, Ug,uji,9) ^ ^ 
dr ' dojg 

The dependence of the critical cluster parameters on the initial molar volume of liquid, 
ui, are shown in Figs. [7H9l for different values of temperature, 9 = 0.17, 0.5, 0.7, 0.8, 0.891 
and 0.92. The positions of the binodal curves are given then by u^'^^^^ = 0.409, 0.45, 0.494, 
0.531, 0.589, 0.62, and = 1.663 ■10^ 90.5, 11.606, 5.634, 3.043, 2.475, the respective 

parts of the spinodal curves are located at u^^^^'^ = 0.452, 0.528, 0.59, 0.634, 0.696, 0.726, 
and wgs^^*) = 13.609, 4.04, 2.584, 2.087, 1.679, 1.547, correspondingly. 

The dependence of the work of formation and radius of the critical cluster on tem- 
perature for the practically significant cases p = —5,0, 1,2,5 and 10 bar is presented in 
Fig. [TOl In Fig. [TTl the dependence of the nucleation rates on temperature for the same 
values of pressure are shown (a value of the pre-exponential factor of Jo = lO^^s^^m^'^ has 
been used for the calculations). One can see, that in such case homogeneous nucleation 
is possible only at very high temperatures, near Tg ~ 1570 K. One can observed as well 
that concerning a 20 cm diameter sphere (region of proton adsorption) and 2 ms time 
span, one can expect 1 or more nucleation event above 1530.5 K. However, heterogeneous 
nucleation may occur also at lower temperatures. 
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5 Conclusions 



In spallation neutron sources, liquid mercury is the subject of big thermal and pressure 
shocks (including negative part), upon adsorbing the proton beam. Increased temperature 
and decreased pressure can cause instable bubbles which can cause cavitation erosion of 
the structural material, shortening the life-time of the equipment and contaminating 
the mercury with tiny steel pieces. Therefore it is crucial to avoid or minimize bubble 
nucleation. While pressure shock can be softened by adding helium micro-bubbles to the 
mercury, there is no way to deal with the thermal shock (i.e. local heating is not possible 
in the middle of the liquid mercury). Therefore our calculation focused on to calculate the 
extent of the temperature increase, the work of critical cluster formation (i.e. the nucleus 
of a macroscopic bubble) and the nucleation rate. It has been shown that after repeated 
proton pulses the temperature can be increased with a few hundred K, but the nucleation 
rate is so low that the possibility of homogeneous nucleation (i.e. bubble formation in 
the pure mercury) is highly improbable, even when the pressure goes below the vapor 
pressure. 
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Figure 1: Isotherms of mercury as described via Eq.([5]) for different values of the reduced 
temperature. 




Figure 2: Comparison of experimental data (according to [131 E]) fo^' the vapor-liquid 
coexistence properties of mercury with the theoretical results (full curve determined via 
Eq. (E])) obtained in this work. 
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Figure 3: Comparison of experimental vapor pressure for mercury (according to [131 [S]) 
with theoretical results (full curve determined via Eq.([5])) obtained in this work. 
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Figure 4: Dependence of the surface tension on temperature, solid curve presents Eq.(|2 
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Figure 6: Gibbs free energy surface for metastable initial state, 6 = 0.92, ui = 0.65. 




Figure 7: Dependence of the critical cluster radius, = Rc/Ra, on the initial molar 
volume of liquid, ui, for different values of temperature, 9 = 0.17, 0.5, 0.7, 0.8, 0.891 and 
0.92. 
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Figure 8: Dependence of the gas molar volume in critical bubble, Wg^c, on the initial molar 
volume of liquid, cji, for different values of temperature, 6 = 0.17, 0.5, 0.7, 0.8, 0.891 and 
0.92. 
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Figure 9: Dependence of the work of critical cluster formation, AGc/kBTQi, on the 
initial molar volume of liquid, uy, for different values of temperature, 9 = 0.17, 0.5, 0.7, 
0.8, 0.891 and 0.92. 
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Figure 10: Dependence of the work of critical cluster formation, I\GcI^bT (a), and of 
the critical cluster radius (b) on temperature for p = —5, 0, 1, 2, 5 and 10 bar. 
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Figure 11: Dependence of the nucleation rate on temperature for p = —5, 0, 1, 2, 5 and 



10 bar. 
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